function [ fx fy ft fxx fyy fxy fxt fyt  ux, uy, vx, vy] = HSGradients(image1, image2, u, v, params)
%     if (nargin < 3)
%         params.spatialDerivBlend = 0.6;
%     end        
    
%     horizKernel = -(1/12) * [ -1, 8, 0 ,-8, 1  ]; % 2nd order central difference scheme    
%     vertKernel = horizKernel';
%     
%     I1x = imfilter(image1, horizKernel, 'symmetric');
%     I1y = imfilter(image1, vertKernel, 'symmetric');    
%     
%     I2x = imfilter(image2, horizKernel, 'symmetric');
%     I2y = imfilter(image2, vertKernel, 'symmetric');        
%     
%     ux = imfilter(u, horizKernel, 'symmetric');
%     uy = imfilter(u, vertKernel, 'symmetric');
%     vx = imfilter(v, horizKernel, 'symmetric');
%     vy = imfilter(v, vertKernel, 'symmetric');
%     
%     ft = image2 - image1;
%     
%     %fx = params.spatialDerivBlend * I2x + (1 - params.spatialDerivBlend) * I1x;
%     %fy = params.spatialDerivBlend * I2y + (1 - params.spatialDerivBlend) * I1y;
%     
%     fx = I2x;
%     fy = I2y;
%     
%     fxx = imfilter(I2x, horizKernel, 'symmetric');
%     fyy = imfilter(I2y, vertKernel, 'symmetric');
%     fxy = imfilter(I2x, vertKernel, 'symmetric');
%     fxt = I2x - I1x;
%     fyt = I2y - I1y;

    [fx, fy, fxx, fyy, fxy] = MyDerivative(image2);
    [f1x, f1y ] = MyDerivative(image1);
    fxt = fx - f1x;
    fyt = fy - f1y;
    ft = image2 - image1;   
    
    [ ux uy ] = MyDerivative(u);
    [ vx vy ] = MyDerivative(v);
    
    fx(isnan(fx)) = 0;
    fy(isnan(fy)) = 0;
    ft(isnan(ft)) = 0;
    fxx(isnan(fxx)) = 0;
    fyy(isnan(fyy)) = 0;
    fxy(isnan(fxy)) = 0;
    fxt(isnan(fxt)) = 0;
    fyt(isnan(fyt)) = 0;    
end

function [ fx fy fxx, fyy, fxy ] = MyDerivative(image)
    [fx, fy, fxx, fyy, fxy] = derivative5(image, 'x', 'y', 'xx', 'yy', 'xy');
    fx(isnan(fx)) = imfilter(image(isnan(fx)), [ -1, 1 ], 'symmetric');
    fy(isnan(fy)) = imfilter(image(isnan(fy)), [ -1 ; 1 ], 'symmetric');    
end